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Abstract 



We present algorithms to compute the Smith Normal Form of ma- 
trices over two families of local rings. The algorithms use the black-box 
model which is suitable for sparse and structured matrices. The algo- 
rithms depend on a number of tools, such as matrix rank computation 
over finite fields, for which the best-known time- and memory-efficient 
algorithms are probabilistic. 

For an n X n matrix A over the ring f[z]/{f^), where is a power 
of an irreducible polynomial / G F[2:] of degree d, our algorithm re- 
quires 0{r]de^n) operations in F, where our black-box is assumed to 
require 0{rf) operations in F to compute a matrix-vector product by 
a vector over F[z]/(/^) (and r/ is assumed greater than nde). The al- 
gorithm only requires additional storage for 0{nde) elements of F. In 
particular, if = 0~{nde), then our algorithm requires only Cr{r?Se') 
operations in F, which is an improvement on known dense methods 
for small d and e. 

For the ring TLjp^TL^ where p is a prime, we give an algorithm which 
is time- and memory-efficient when the number of nontrivial invariant 
factors is small. We describe a method for dimension reduction while 
preserving the invariant factors. The time complexity is essentially 
linear in /j^nrelogp, where is the number of operations in Z/pZ to 
evaluate the black-box (assumed greater than n) and r is the total 
number of non-zero invariant factors. To avoid the practical cost of 
conditioning, we give a Monte Carlo certificate, which at low cost, 
provides either a high probability of success or a proof of failure. The 
quest for a time- and memory-efficient solution without restrictions 
on the number of nontrivial invariant factors remains open. We offer 
a conjecture which may contribute toward that end. 

Category: G.4. Mathematical SoftwareAlgorithm Design and Analysis 
Category: 1.1.4. Symbolic and Algebraic ManipulationApplications 
Terms: Algorithms, Complexity, Performance. 

Keywords: Local Principal Ideal Ring, Sparse Matrix, Polynomial Matrix, 
Integer Matrix, Smith Form, Black Box. 
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1 Introduction 



We consider the problem of computing the Smith Normal Form (SNF) of 
sparse matrices over (commutative) local principal ideal rings (PIRs). The 
Smith form is a diagonalization of matrices which has many applications 
in diophantine analysis (Chou and Collins, 1982), integer programming (Hu, 
1969), combinatorics (Wallis ct al., 1972), determining the structure of Abelian 
groups (Newman, 1972) and class groups (Hafner and McCurlcy, 1989), com- 
puting Simplicial Homology (Dumas et al., 2003), in system theory (Kailath, 
1980; McMillan, 1952), and in the study of symplectic spaces (Chandler et al., 
2010). 

The original work of Smith (1861) proved existence and uniqueness of the 
SNF for integer matrices. The generalization to PIRs is due to Kaplansky 
(1949). 

The problem of computing the Smith form of a sparse matrix over a 
principal ideal ring presents several challenges. One approach is to simply 
compute the SNF over the global ring (i.e., F[z] or Z) and then reduce the re- 
sult modulo the power of the prime ideal. The algorithm of Gicsbrccht (2001) 
for SNF of a sparse matrix over Z could be used, but the ultimate time re- 
quirement is essentially cubic (although space requirements are lower). An 
asymptotically faster algorithm along similar lines, but which requires con- 
siderably more space, is presented in (Eberly et al., 2007). The best known 
algorithm for dense matrices over F[z] by (Storjohann, 2000, Prop 7.16) re- 
quires time essentially equal to matrix multiplication. However, it is not 
sensitive to sparsity. 

On the other hand, computations over F[z] and Z suffer from coefficient 
growth which is not clearly necessary in a PIR. For example, over Z/p^Z 
where p is a prime, one might hope to perform all computations modulo 
p^, and not with integers larger than p^. Storjohann (2003) provides a fast 
algorithm using elimination in a PIR, but it is not sensitive to sparsity and 
requires time proportional to matrix multiplication. Wilkening and Yu (2011) 
demonstrates an algorithm for dense polynomial matrices over local rings, but 
offers no complexity analysis. Dumas et al. (2001) give black-box algorithms 
over Z and locally at a prime, which however do not have a benefit when 
only a few invariant factors are nontrivial. 

When dealing with sparse matrices we would like to preserve the spar- 
sity of the input matrix, and introduce no fill-in. Thus we pursue black-box 
algorithms in the sense that the input matrix is only used for matrix-vector 
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products. The complexity of black-box algorithms is thus expressed in terms 
of the number of matrix-vector products used. Space requirements are kept 
to the storage of a few vectors. There has been great success in applying 
black-box methods over finite and arbitrary fields, starting with Wiedemann 
(1986), where the cost of many linear algebra problems has been reduced to 
computing a linear number of matrix-vector products. Our ultimate goal is 
then to add local Smith form to that list. 

Specifically we will consider the local Artinian principal rings (also known 
as special principal rings) Z/p'^Z for a prime p and positive exponent e, and 
F[z]//'^F[2;], for irreducible / G f[z]. Let L be a local Artinian principal ideal 
ring with a maximal prime ideal pL. For any matrix A G L"^", there exist 
unimodular matrices U,V & \_nxn ^ diagonal matrix S G L"^" such that 
A = USV, where 



Definition 1. 5* is called the Smith form of A, and the diagonal elements 
are called the invariant factors of A. 

Our goal throughout this paper, is to compute the multiplicities of the 
Smith form invariants, i.e, {tq, ri, . . . , re_i} for given black-box matrices, and 
in particular, sparse matrices. 

Our Contribution. For matrices over F[z]/{f'^), we present an algorithm 
which relies on computing ranks of related black-box matrices over F, and 
give a complete complexity analysis. The key idea of our algorithm is a linear 
representation of polynomials in the ring F[z]/{f'^) as matrices over F, and 
using rank computations over F to discover the multiplicities of the Smith 
invariants. This reduction allows us to take advantage of the well-studied 
efficient algorithms for computing ranks of sparse matrices over fields rather 
than rings. The cost of our algorithm is 0{'r]de'^n) operations in F, where each 
black-box evaluation costs t] operations. Our approach takes a path similar 
to the linearization of Kaltofen et al. (1990) for matrices over F[z]. This 
approach is also explored for dense matrices over local rings by Wilkcning 
and Yu (2011). Dumas ct al. (2009) used rank computations to discover 
multiplicities of characteristic polynomial factors for black-box matrices over 
fields. 

The linearization idea, however, would not be applicable over the integers 
since there are no appropriate linear representations from Z/p'^Z to Z"^". 



S = diag(l, . . . , l,p, . . . ,p, . . . . . . ,/-\0, . . . , 0) 

— — — — ■' ^ ■' 




(1) 
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Hence, it is necessary to develop different methods for the integer case. A 
useful approach to Smith form computation is to begin by determining which 
primes occur in the invariant factors and then compute the form locally 
at those primes. This has been done in several recent algorithms (Dumas 
et al., 2003; Saunders and Wan, 2004). A fully memory efficient, black-box 
algorithm for sparse and structured matrices has not been given, however, 
for lack of an efficient black-box algorithm for the Smith form locally at a 
prime. 

Toward that end we give a black-box algorithm over Z/p^Z, whose cost 
is essentially dominated by unek for an n x n sparse matrix with /i nonzero 
entries, e being the largest exponent of the prime in the Smith form, and k = 
YliZi being the number of nontrivial invariant factors (ro is not included). 
It is to be expected that the cost depends on both fin and e. The dependence 
on k, although unfortunate, is not completely unlikely. It is natural that it 
is easier to find Smith form for matrices with fewer number of non-trivial 
factors. In addition, both e and k are small in many cases of interest. Some 
applications with this property are discussed at the beginning of Section 3. 
The key idea of this algorithm is to apply a reduction in dimension to dispose 
the zero invariant factors, compute a nullspace basis of the reduced matrix 
to dispose the ones, and then determine the nontrivial invariant factors by 
dense elimination methods. This idea is applicable to the polynomial case as 
well. However, it is not interesting, since the complexity of this method has 
an extra factor of k over the rank-based method presented in Section 2. 

Notation. Throughout this paper F denotes a field and L denotes a local 
ring. We use Zp to denote Z/pZ and Z^e to denote Z/p'^Z. In the complexity 
analysis, we count the algebraic complexity in the base field, i.e. we assume 
that all operations in the base field have a unit cost. We use "soft-Oh" to 
hide the logarithmic factors. We say that / G 0~{g) if there exists a constant 
c such that / G 0{g log'^ g). We use M(n) to denote the number of operations 
in the base field required to multiply two polynomials of degree at most n. 
Finally, we use 0{n'^) to denote the matrix multiplication exponent, e.g., 
u < 2.372 using Coppersmith and Winograd (1990). 

Roadmap. The rest of this paper is organized as follows. Section 2 con- 
tains the algorithm and analysis for Smith form over F[z]/(/^). In Section 3, 
our nullspace algorithm for Smith form over Zpe is given after some discussion 
of applications, a development of preconditioners for the problem. A Monte 
Carlo verification method is given that can be of use with small primes. A 
conjecture is also given, that may shed some light on the problem when there 
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are many nontrivial invariants. Finally, Section 4 is a brief summary. 



2 Smith form over F[2;]/(/^) 

In this section we present an algorithm to compute the local Smith form of 
a sparse polynomial matrix. Throughout this section, let F be a field, L = 
where / G F[z] is irreducible of degree d, and e G Z>i. The ideals 
in this ring are of the form /*L for < i < e and the RHS of equation (1) 
becomes 

diag(l^_^, /^_^, . . . , r-\...,r-\ 0, . . . , 0) (2) 

ro ri re-1 

Our goal is efficiently compute the multiplicities: {ro, ri, re_i}. The 
approach is to embed the ring L"^" in the ring p"'^^^'^'^'^ and reduce the 
computation to finding ranks of matrices in the base field, F, where known 
fast black-box algorithms can be used. 



2.1 Embedding of L^^" in F"^^^"^^ 

In this section, we describe the classical embedding of L into F'^^^'^^, and 
how properties of matrices over L are revealed by their images over F. First, 
define the map '■ ^ ^ fdexde follows. Suppose = + aiz + ■ ■ ■ + 
a(ie-iz'^^~^ + z'^'^, with a companion matrix 



Ct 



/O 

1 ••. 
Vo ■■■ 



-ao \ 

-ai 

-O.de-1/ 



Define ^e{z) = Cfe, and ^Pe{z^) = ^eizY- By linearity, extend to all 
elements oi g = + giz + ■ ■ ■ + gde-iz'^'^~^ G L such that ipe{g) G F*^^^*^^: 



^eig) = giCfe) = gol + giCfe + g2C'}e H h gde~iC[ 



de-l 



It is straightforward to verify that (fe is a ring isomorphism between L and 
F[C,p]. 



Lemma 2. rank(y9e(/*)) = d{e — i) for < i < e. 
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Proof. Since f{Cfe) acts as multiplication by / mod f^, it has null vectors 
which are images of polynomials in /^~^L. This is a vector space of dimension 
d, and hence rank(/(Cfe)) = de — d. Also, f{Cfe) has minimal polynomial 
x^, whence 

/Od Id ■■■ Od\ 

' Id 

\0d OdJ 

where Id, Od € F*^^*^ are identity and zero matrices respectively. The rank 
ifeif^) = feifY is now evident from the structure of the nilpotent Nf. □ 

We extend the map <y9e to n x n matrices over L. For every A G L"^", 
ipe{A) is a nde x nde matrix over F, where every entry aij of A is replaced 
by the de x de block ipeidij)- Applying ip^ to (2), we get 

iPe{S) = diag(v3e(l), • • • , V5e(l), ^df), • • • , ^eif), ■ ■ ■ , 
V ' ^ V ' 

ro ri 

y,e(r-'),...,¥^e(r-'),0,...,0)GF"'^-"'^^ (3) 
V ' 

re-1 

Lemma 3. If U G L"^" is invertible, then ipe{U) G F"'^'^^"''^ is invertible. 

Proof. If U is invertible, then there exists a G L"^" such that t/VT = /, 
and (fe{U)(fe{W) = (pe{I)- But (fe{I) = Inde, SO (fe{U) has iuvcrsc V5e(W^)- □ 

We now establish the property relating multiplicities in the invariant fac- 
tors of A G L"^" to the rank of v^eC-)- 

Theorem 4. Let A G L"^" have Smith form in (2), then Taiik{(fe{A)) = 
devQ + d{e — l)ri + ■ ■ ■ + dr^-i. 

Proof. There exist unimodular matrices U,V & \_nxn g^j^}^ that UAV = S. 
By isomorphism, ipe{U)(fe{A)ipe{V) = ipe{S). By Lemma 3, ipe{U),(pe{V) 
are invertible and thus Ya.nk{ipe{A)) = rank(<y9e(S')). In (3), (pe{S) is a block 
diagonal matrix, so 

e-l 

rank(v9e(5')) = ^ rank(v9e (/')). 

i=0 

The proof then follows from Lemma 2. □ 
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A consequence of Theorem 4 is that, for 1 < £ < e, we have rank(y9£(A mod 
/^)) = diro + d{i — l)ri + - ■ ■ + dre-i. For example, lank^ipi^A mod /)) = drQ. 
In general, we have the following corollary. The proof is left to the reader. 

Corollary 5. Let p£_i denote rank(y9e(^ mod /^)), 1 < i < e. Then 



( d 

2d 

\ed 





d 



2d 



0\ 



dj 



/ro\ 

ri 

\Te-lj 



Pl 
\Pe-lJ 



(4) 



This system may be solved in linear time. Let ai = pi and 0"^ = pi — Pi~i, 
for 1 < i < e. Then the cTj are the prefix sums of the ri, so ri = cxi and 
'"i = o"i — cTj-i, for 1 < i < e. 

Next we consider how to efficiently compute {po, pi, . . . , Pe-i} for a given 
black-box matrix. 



2.2 Black- box for the embedding 

Given a black-box for A G L"^", over L = f[z]/{f'^), we can easily construct 
a black-box for V5£(v4 mod /^), for all £ < e, at not much higher cost. 

First, we define the black-box model cost model, and then show how to 
perform black-box computations under transformations efficiently. 

Definition 6. Let A G L"^" be a sparse matrix over L = ^[z]/{f'^), where 
f G F[z] is monic and irreducible of degree d. The black-box for A is a 
mapping L" — )■ L" such that for all f G L", Av G L" can be computed with r) 
operations in F. We assume throughout that t] > nde. 

Lemma 7. Suppose we are given a black-box for A G L"^", where L = 
^lA/if^) above. Let £ G {1, . . . , e} and v G F'^^" with unique pre-image v G 
f[z]/{f^). Then we can compute (pe{A mod f^)v G F'^^" with 0{r] + n M{de)) 
operations in F. 

Proof. Assume that v G F*^^" is labelled as: 

V = {Vifi, . . . ,Vi^de-l,V2,0, ■ ■ ■,V2,de-l, ■ ■ ■ ,Vn,0, ■ ■ ■,Vn,de-l)- 

Construct the vector v = {vi, . . . G L", where Vi = J2o<j<de^i,j^'' ^ '"W- 
Now, compute w = Av mod G using rj operations for the black-box 
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evaluation plus 0{nM{de)) operations in F. Let w = {wi, . . . ,Wn)- Assume 
= J2o<j<di^i,jZ^ e F[z]. Then 

□ 

Our algorithm for computing the Smith form of a matrix A G L"^" given 
by a black-box is now straightforward. Using Theorem 4 and Lemma 7 we 
reduce the computation of p^'s in (4) to computing ranks of matrices over the 
ground field F, which can be accomplished using existing fast and memory- 
efficient black-box algorithms over fields, e.g. Wiedemann's algorithm. 

Algorithms for computing the rank of a black-box matrix over a field are 
developed by Wiedemann (1986), and refined in subsequent work of Kaltofen 
and Saunders (1991), Ebcrly (2004), and others. If the input matrix is in 
pnxn g^^^ ^YiQ black-box evaluation requires r] operations in F, then the rank 
algorithms require 0~{nri) operations in F. They are probabilistic, and return 
the correct rank with controUably high probability on any input. We will 
assume that an appropriate choice of black-box rank method is made, and 
note that there is considerable difference in their effectiveness in practice and 
over various ground fields. 

Algorithm 1. Smith invariants in L = F[2;]/(/'^), where f G F[z] is irre- 
ducible of degree d. 
Input: Black-box for A G L"^". 

Output: To, . . . ,re_i such that ri is the multiplicity of /* in the Smith Form 
of A, and the multiplicity of is n — J2i ^i- 

1. For all i E {1, . . . , e}, invoke a black-box rank algorithm on the black- 
box for ipe{A mod /^) : F'^^" — > F'^^". Let pg_i = Tank{(fi{A mod /^)). 

2. Solve (4) for tq, . . . , r^-i- 

3. Return tq, . . . , re_i. 

Theorem 8. Algorithm 1 is correct, and requires 0~{r)de'^n) operations in F. 
The space requirement of the algorithm is 0{den) elements in F. 

Proof. The correctness follows from the results and discussion in this section. 
We analyze the time and space complexity of step (1), which dominates. This 
step requires O^de^n) black-box evaluations, and storage for 0{nde) elements 
in F. □ 
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3 Smith form over Z 



In our experience, most Smith normal forms of integer matrices in practice 
involve relatively few non-trivial factors, i.e., most of the invariant factors 
are I's or O's. The algorithm of this section addresses that situation. 

For example, in some work on computing homology of simplicial com- 
plexes, Dumas ct al. (2000, 2001); Babson et al. (1999); Bjorncr and Welker 
(1999), large boundary matrices arose. One of the most challenging Smith 
forms to compute at the time was a 135135 by 270270 matrix which turned 
out to have 133991 ones, 220 3's, and 924 zeroes as the invariants. Other 
examples in that study were also large but with even fewer nontrivial in- 
variants. Most often in homology computation, the number of I's (the Betti 
number) greatly exceeds the number of non-trivial entries, it seems. 

For another example, recently in the study of symplectic 3 spaces. Smith 
form computations have been desired of some rather large matrices Chandler 
et al. (2010); Chandler (2011). For these matrices it is conjectured that only 
one prime will appear in the invariant factors (other than the largest) and 
indeed only the local Smith form at that prime is of interest. Furthermore, 
the conjectured structure predicts only a few nontrivial invariant factors. We 
denote these examples as W3-g, where g is a prime power and the Smith form 
modulo q is desired. W3-g is a {0,l}-matrix of size approximately x q^ 
with about q"^ nonzero entries. The current challenge is to compute Smith 
form of W3-64 and W3-81. The algorithm described here is designed to 
handle this case. 

For a prime p and an exponent e G let ip be the natural projection 
Zpe — )• Zp, which extends naturally to ip : Z^i^" — )■ Z^^" by element-wise 
mapping. Note that x G Z^e is a unit if and only if '{'{x) ^ 0. Likewise, 
A G Zpe^" is unimodular if and only if (p{A) is unimodular. 

3.1 Nullspace method 

Let us introduce the approach by way of a sketched example. Suppose A is 
a matrix over Zps. Further suppose A is 100 by 100 and 

A~ diag(l,l,...,l,p,p,p,p^p^O,0,...,0), 

with 45 ones and 50 zeroes. This approximates on a small scale the pattern of 
invariants we expect to see on W3-g. First, a reduction in dimension allows 
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us to reduce A to an i x i matrix p{A) having the same nonzero invariants, 
where i is the (max) rank, or shghtly larger. We illustrate with i = 52: 

p{A) ~ 5 = diag(l, 1, . . . , 0, 0). 

Over Zp, the nullspace basis N' of (p{p{A)) {N' has 7 columns) is equiv- 
alent to that of S. Let E' be the last 7 columns of the 52 x 52 identity 
matrix. Let E and N be arbitrary embeddings of E' and N' in Zp^^, such 
that (p{E) = E',ip{N) = N'. Thus p{A)N and 5*^ are multiples of p and 

piA)N r^SE = diag(p,p,p,p^p^O,0). 

In summary, the algorithm is to apply a reduction in dimension to dispose 
of zeroes, compute nullspace basis to dispose of ones, and determine the 
nontrivial invariants by computing Smith form of AN using dense methods. 
AN is an n X A; matrix, where k is the number of nontrivial invariants (or 
slightly larger - i - rank mod 2). 

Reduction in dimension is a frequent tool and has been used for Smith 
form, for example, in Dumas et al. (2001). But then their computation 
proceeds without disposing of the unit invariant factors. Thus the time 
complexities below, otherwise similar to theirs, differ in that we replace a 
rank factor £ by the number of nontrivial invariants, k. 

3.2 Probabilistic dimension reduction 

Let A G Zpc^", for which we have a fast black-box. Let A have a Smith form 
diag(si, . . . , Sr, 0, . . . , 0) G Z^e^". Our goal in this subsection is to construct 
p{A). That is, given such an A and a £ G {1, . . . , n}, to construct a black-box 
of similar cost for a matrix B G Z^?^ which has Smith form diag(si, . . . ,Si), 
i.e., with the initial invariant factors of A. 

Notationally, for integers n and k < n, let be the set of fc-tuples of 
distinct elements (in increasing order) of {1, ... , n}. For a matrix B G L"^", 
and a,T E C^, define B(^^ as the (cr, r) minor of B, i.e., the determinant of 
the k X k submatrix of B with rows from cr and columns from r. We use 
script letters, e.g. 2),T, to denote matrices with indeterminate entries. 

We use techniques similar to that derived in Giesbrecht (2001) with scaled 
Toeplitz matrix conditioners. For indeterminates A = {vi,Wi,yi}, let Di = 
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diag(wi, . . . , Vn), = diag(wi, . . . , Wn), and T be a generic Toeplitz matrix: 

/ Vn Vn+i ■■■ yi \ 



■ Vn ■■ ■ 

yn-2 ' ■ • Vn+l 

\y2n~l l/2n-2 ■■■ Vn / 



(5) 



Lemma 9. Let ^ = 2)iT2)2 in- the indeterminates A, as in (5). Let k e 
{1, . . . ,n} and a = (ai, . . . , ak) , T = {ti, . . . , r^) G 

T(^) G Z[A] /ias contend 1; 

25 (^) = ■ ■ ■ i^^feW^i ■ ■ ■ . 

Proof. Part (i) is from (Giesbrecht, 2001, Lemma 1.3) and part (ii) follows 
easily from the Cauchy-Binet formula. □ 

Note that (|^) uniquely identifies which minor of ?B was selected. 

Lemma 10. Let A G Z"^", and 25i, be n x n matrices of distinct in- 
determinates from a set A, of the form (5). Then 21 = ^iA^2 is such 
that for 1 < k < n, the content of i/jk = 2t(J"'^) G Z[A] equals Ak, the kth 
determinantal divisor of A. 

Proof. By the Cauchy-Binet formula we have 



i...k 



Thus 2l(|^' "1^) is a sum of polynomials of content 1, with distinct indetermi- 
nates, one for each kxk minor of A, times the value of that minor. Hence it 
must have content equal to the GCD of all /c x /c minors of A, which is equal 
to the fcth determinantal divisor. □ 

Theorem 11. Let A G Z"^", p > Qri^^ a prime, and ^ > 2. Let -81,-82 G 
^nxra f)Q formed by a random assignment of variables in Q3i,232 in (5) respec- 
tively, where choices are made uniformly from L = {0, . . . ,6n^,^ — 1}, and 
A = B1AB2. Then with probability at least 1 — 1/^, for all 1 < A; < ri, the 
order of p in A^, the kth determinantal divisor of A equals the order of p in 
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Proof. Let ipk be as in Lemma 10, which has content equal to the kth determi- 
nantal divisor of A. Observe from our construction that deg ipk Qk < 6n 
(the total degree of a kx k minor of an indeterminate Toeplitz is < k). Thus, 
with values selected as described, by the Schwartz-Zippel Lemma (Zippcl, 
1979; Schwartz, 1980), we have that ipk/^k is a polynomial in the entries of 
matrices Bi , B2 and 

prob{(^/'A:/Afc) ^ Omodp} > 1 



This is exactly the probability that the order of p in equals the order of 
p in the leading k x k minor of A. The probability that this happens for all 
k, from 1 < A; < n is at least (1 - l/«)))" > 1 - 1/^. □ 

Corollary 12. Let p > Qn'^C, be prime, for a ^ > 1, and e > 1, and suppose 
A e Zpi"" has (local) Smith form diag(si, . . . , s„) G Zpi"". Let Bi, B2 € 
Zpc^" be formed by a random assignments of variables in 53i,532 G Z"^" in 
(5) respectively, where choices are made uniformly from L = {0, . . . , 6n^^ — 
1} modp^ Let A = B1AB2 G Z^i"", and for 1 < k < n let Ak be the 
leading kx k submatrix of A. Then with probability at least 1 — 1/^, for all 
k G {1, . . . , n], the local Smith form of A^ is diag(si, . . . , s^) G Zpe^^. 

Proof. The Smith form of A equals the Smith form of any A G Z"^" with 
A = A mod p*^, reduced modulo p'^ (the only non-units will be powers of p 
after the reduction). Thus, Theorem 11 implies that the order of p in the 
kth determinantal divisor of A equals the order of p in the leading k x k 
minor of A, for all k, with probability at least 1 — 1/.^. This implies that 
y4fc = Ak mod p"^ will have Smith form (si, . . . , Sk) for all 1 < A; < n where 
Ak is the leading k x k minor of A, since A^ = si ■ ■ ■ Sk for 1 < k < n. □ 

Computationally, if we know that rank of A is at most m, then we can 
work with truncated random scaled Toeplitz matrices Bi G Z^^" and B2 G 
^nxm_ r^Yien Corollary 12 implies that A = B1AB2 G Z^e"""" has the same 
non-zero invariant factors as A. 

Working with Small Primes. The conditions for Corollary 12 require 
p > 6n^^. For smaller primes the algorithm may well work, but appears 
much more difficult to prove. The following method may be used to remedy 
this. 



13 



The approach is to replace Z by a subring of the ring of algebraic integers 
in a number field of degree 1] = \logp{6n'^^)~\ over Q, in which p is inert. 
Specifically, let T G Z[y] have degree at least r] be such that F mod p is 
irreducible in Let 7 e C be a root of T. Then ^[7] is such that 

Z[7]/(p^) is a local ring which contains Z^e, and such that the residue class 
field Z[7]/(p) contains more than Gn^^ elements. We call Z[7]/(p^) the Galois 
ring with p^ elements, and denote it by GR{p'^,ri) (see McDonald (1974)). 
Like Zpe, GR(j9*^,?7) is a local principal ideal ring with maximal prime ideal 
generated by p. 

Analogues of Theorem 11 (over ^[7]) and Corollary 12 (over GR(p^,?7)) 
can be proven similarly. We state the latter formally, but leave the proofs to 
the reader. 

Corollary 13. Let p be prime, e > 1, > 1 and rj = [logp(6n^^)] . Let 
GR{p'^,ri) = Z[y]/(r) for T G Z[y] of degree rj which is irreducible modulo 
p. Suppose A G Zpe^ has (local) Smith form diag(si, . . . , s„) G Zpi^". Let 
-81,-82 G GR(p, T])"^" be formed by random assignments of indeterminates 
in !Bi,232 in (5) respectively, where choices are made uniformly from L = 
{Zo<i<r,(^if ■ e {0,...,p- l}}modp^ LetA = B^AB2 G GR(p,r/)"X", 
and for 1 < k < n let Ak be the leading k x k submatrix of A. Then with 
probability at least 1 — 1/,^, for all A; G {1, ... , n}, the local Smith form of Ak 
IS diag(si,...,Sfc) G Zle^. 

3.3 Probabilistic validation of dimension reduction 

This section might be titled, "Escaping the tyranny of Schwartz-Zippel." 
It is frequently the case, as is exemplified in the previous section, that an 
argument for a favourable probability based on the Schwartz-Zippel Lemma 
requires an inconveniently large set for random assignments. Experience in 
practice demonstrates that far smaller sets suffice virtually always. In this 
section we give a method to have a provably low probability of failure while 
making no assumptions about the basic preconditioner. 

For a matrix A let Sk{A) denote the k-ih. invariant factor of A (in the order 
in which Si{A) \ Sj+i(A), for 1 < i < n). Let Sk{A) denote the leading k x k 
submatrix of the Smith form of A, Sk{A) = dia.g{si{A), S2{A), . . . , Sk{A)). 
Let [A, B] denote the side by side join of two conformable matrices. 

Theorem 14. Let A G L'"^", Q G L"^'', and y G L", with k < min(m, tt.) 
and the entries of y being uniform random variables over L. Then AQ has 
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the first k invariant factors of A with high probability if AQ and A[Q, y] 
have the same first k invariant factors. To be precise, we have the following 
conditional probability: 

prob(^,(A) = Sk{AQ) I Sk{AQ) = Sk{[AQ,Ay])) > 1 - 1/p. 

Proof. Note that Sk{A) = Sk{B) if and only if the i-th detcrminantal divisors 
are equal: Di{A) = Di{B),i G l,...,k, where Di{A) denotes the GCD of 
all /l(^), cr, r G C". We will call any i x i minor equal to Di{A) (up to unit 
multiple) a witness for the determinantal divisor Di{A). Note that the GCD 
of a set of elements of L is a member of the set (up to unit multiple) , so every 
determinantal divisor has at least one witness. For example, for a Smith 
form S, the j-th determinantal divisor, if not a unit, has exactly one witness, 
Dj{S) = ^CgJ), where we define l by = (1, 2, . . . 

Without loss of generality, suppose for notational symplicity that k < 
n < m. Let A — USV be the Smith form factorization of A with S — 
diag(si, . . . , Sn) and f/, V unimodular. We are concerned with the invariants 
of A,AQ, and Multiplication by unimodular does not affect 

invariants, so let us consider SV,SVQ, and S[VQ,Vy]. Because V is uni- 
modular, Vy is a uniform random variable if y is. Also, we have made no 
conditions on Q, so we simplify notation, substituting Q for VQ and y for 
Vy. In other words, without loss of generality we may assume ^ = ^S" is in 
Smith form. 

We will show the contrapositive of our proposition. Suppose j is the 
first index at which Sj{A) ^ Sj{AQ). Then ^(^gj) 7^ A(;g)Q('(j)), for all 

(7 G Cj. Otherwise such a minor would witness equality of Sj{A) and Sj{AQ). 

It follows that p I QCi^^). 

Because j is the first such case, there must be an (j — 1) x (j — 1) minor 
Q{'^^''^^^) which is a unit (not divisible by p), where i(j — 1), r G Let r' 

denote rU{A;+l}, a column index set for [Q, y] which includes the last column. 
Then the expansion of the i by r' minor of A[Q, y] has as coefficient of yj the 
term A{^^^^}^)Q('^^^^^) , which is a unit. Thus, for each setting of yi, . . . ,yj-i, 

there is at most one value modulo p or yj for which p \ [Q,y]{''^f') ■ For all 
other values, this minor witnesses that Dj{AQ) ^ £)j(A[(5,y]. Thus when 
S'fc(A) 7^ S}i{AQ) there is at most a 1/p chance that 5'fc(A[(5, y]) agrees with 
Sk{AQ). □ 

The following corollary allows us to verify or disprove the success of pre- 
conditioners such as those used in dimension reduction. This works when the 
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theory of the preconditioner's probabihty of success is invahdated because of 
an insufficiently large set from which random values are chosen. Thus in 
practice for a small prime, one can skip the domain extension described at 
the end of Section 3.2. Indeed, no theory of the preconditioner is required 
at all. One can try to "get lucky" and skip preconditoners altogether. For 
instance, apply the corollary where PAQ selects the leading k x k submatrix 
of A. If the method validates, then the Smith form was found economically, 
otherwise try a more thorough preconditioning. 

Corollary 15 (projection verification). Let matrix A G Z"^" be given and 
heuristic preconditioners P E T}^^ and Q E TI^^^ . Choose Ri E 77^^'^, R2 E 

at random. Let B = ( ^^^^ ^^!^rl] e Z^+'=^'=+^ If S J PAQ) = 

\K2AQ K2AK1 J 

Sk{B) then these are the first k invariant factors of A with probability greater 
than 1 — 2/p'^. 

f PA\ 

Proof. Let C = ij^j^j so that B = [CQ,CRi]. Note that the condition 

Sk{PAQ) = Sk{B) implies that all minors of B containing PAQ have those 
invariants. In particular B and its left side CQ have the same invariants. So 
Theorem 14 applies c times for each of the columns of -Ri. Because these are 
independent random vectors the probability that SkiCQ) = Sk{[CQ,CRi]) 
when Sk{CQ) 7^ Sk{C) is at most 1/p'^. Then apply Theorem 14 to A, PA, C 
on the left c times for the c rows of R2- Again, the probability of an unfor- 
tunate equality is at most l/p^, and otherwise Sk{PA) = Sk{A) is validated. 
Thus the overall probability of success is at least (1 — 1/p'^)'^ > 1 — 2/p'^. □ 

The idea to use some random dense rows in preconditioning is widespread. 
It was used already in (Wiedemann, 1986, proof of theorem 1), and in a 
sense it is the basis of block iterative methods. The idea to obtain a good 
probability (especially for small primes) by solving a sparse problem twice, 
first without the few random dense rows and/or columns then with them, was 
used in Saunders and Yousc (2009). For integer lattices, there are also some 
similarities to the additive preconditioners of Eberly et al. (2000), and the 
lattice compression of Chen and Storjohann (2005), especially in the analysis 
for small primes dividing invariant factors. 

The parameter c can be adjusted to get the desired degree of certainty. 
For example, c = 21 ensures probability of failure less than one in a million, 
since 2/p'^ < 1/2^^ < 10~® in that case. Also, one can pad with more random 
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rows and columns to improve weak preconditioners. Thus if c = 30 is used 
and the matrix with the first 9 of those columns and rows has the same 
Smith form as the matrix with all 30 random rows and 30 random columns 
adjoined, we have computed the Smith form with error expected less than 
once in a million trials. In effect, we have corrected for some weakness in the 
heuristic preconditioners with 9 extra rows and columns and then verified 
with 21 more. 

3.4 Algorithm for the Smith Normal Form 

After reducing the dimension to a value at or near the number of nonzero 
invariant factors, the following algorithm is applied. Recall that is the 
natural projection Z^e Zp. 

Algorithm 2. Smithpe-nuUspace 

Input: a black-box for B e 11^'^ , and a bound i for the number of nonzero 

invariant factors 

Output: S, the Smith form of B. 

0. Set A to the dimension reduction of B to i x i. 

1. Let ro = iai\k{ip{A)) over Zp. The nullity of ip{A) is then k = i — rQ. 

2. Compute N' e Z^?*^, a lifting to Zpe of a right nullspace basis of ip{A) 
over Zp. 

3. Let N = AN' G Z^?^. This involves k matrix vector products with A. 
Note that N is divisible by p. 

4- Compute the Smith normal form of N over Xpe by Gaussian elimina- 
tion: 

diag(p^_^, P^.^.,/ , . . . , f-^,..^.,p^-\ , 0^^^^^ . 

ri r-2 re-i 

5. Return tq, . . . , r^-i- 

We will analyze the algorithm holding e and p constant. Considering 
them as parameters would introduce a factor of O~(elog(p)). Let the cost of 
matrix- vector product by B is 0{fi). Since we arc holding e and p constant, 
this is the same for application to vectors in Z^' and in Z^e. 

Step 0: Toeplitz matrices may be applied to vectors via polynomial mul- 
tiplication, so the cost of the black-box for A is 0{M{n) -\- jj). M{n) is 0~{n) 
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and we will assume /i > n, so the black-box cost of matrix- vector product by 
A is 

Step 1: The rank over Zp can be done by a black-box method in 0~{{ifi) log(.^)) 
to achieve probability of error less than 1/^ (Wiedemann, 1986). Memory 
requirement is 0{1) vectors in Z^. 

Step 2: Let k = i — tq denote the nullity of A modulo p. By black-box 
methods, k random samples of the nuUspace will yield a nullspace basis N'. 
Oversampling can be done and column echelon form computation used to 
reduce to a basis of k columns if need be. The cost is 0~{k{ifi)). Space is 
0{ki). For instance see Chen et al. (2002). 

Step 3: The cost of applying A to N' is 0~{k^). 

Step 4: Any nullspace for 5* over Zp is of the form EW', where E is 
the last i — ro columns of the identity and W is k x k unimodular. Then 
= AN' = USVN = USEW modulo p, for some unimodular W. This lifts 
to a factorization AN = USEW modulo p*^ with U, W unimodular. Thus 
AN has the Smith form SE. The local Smith form of this dense matrix 
can be computed by elimination. An algorithm running in (T{ik'^~^) is in 
Storjohann (2000). 

Theorem 16. Algorithm 2 is a correct Monte Carlo algorithm for computing 
Smith normal form overZpe. The time complexity is 0~{ik{k'^~'^ + fi)) , where 
k is the number of nontrivial (neither nor 1) invariant factors, i is the 
reduced dimension (which can be the rank), and /i is the cost of matrix vector 
product by B. Note that the time complexity is 0~{ik^) under the very modest 
assumption that k'^~'^ < ^. The memory requirement is 0{kt). 

3.5 A conjecture about p-adic carries 

While attempting a general approach for the Zpc case without using dense 
elimination and while using only rank computations over the field Zp, we 
discovered an interesting pattern about ranks of matrices over Zpe which 
could be of independent interest. To put the discovered conjecture in a 
proper context, we first introduce the attempt which lead to discovering it, 
then we present the conjecture which discourages this approach. 

Consider an element a G Zpe written in terms of its unique p-adic expan- 
sion as a = ao + «iP + • • • + ae-iP*^"^ where G Zp for < i < e — 1. This 
representation extends naturally to matrices over Zpe, i.e., A = Aq + pAi + 
. . .-\-p'^~^ Ae-ii where Ai G Z^^". For clarity of presentation and limited space. 
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we only confine the discussion to the case e = 2. Expanding A = USV, we 
get (Ao+pAi) = {Uo+pUi){So+pSi){Vo+pVi), where Aq = UoSqVo mod p 
and 



A 



1 



U.SoVo + UoSoVi + UoS.V, 



(6) 



where rank(S'o) = tq, rank(S'i) = ri, which are the multiplicities of I's and 
p's in the Smith form diagonal, respectively. Furthermore, with appropriate 
preconditioning*, rank(Ao) is proportional to tq, and rank(Ai) is proportional 
to tq + 2ri. Hence, this formulation leads to a belief that we can easily 
isolate Aq,Ai, compute their ranks over Zp, and discover multiplicities of 
the invariant factors. However, by closer inspection, equation 6 is in fact: 



where the extra term, {UoSoVo quo p) , is introduced by the fact that opera- 
tions in Zpe exhibit carries. As a simple example 5-5 over Z34 is (2+p)(2+p) = 
l + 2p + 2p^, where 2p^ is a carry term. These carries contribute to the over- 
all ranks of matrix expressions, and present a challenge in computing the 
Smith form. We hoped to reasonably bound the ranks of these carries, so a 
working algorithm could be developed which reduces Smith form computa- 
tion to efficient rank computations over Zp. This approach would be superior 
to algorithms presented in previous section and literature since it preserves 
sparsity. It is worth noting that in the polynomial case, this approach yields 
a working algorithm. The following conjecture illustrates why the p-adic case 
is more difficult to resolve in this way. 

Conjecture 17. Assume p is a prime, U,S,V G Zp^" such that U, V are 
invertible, S = diag(l, . . . , 1, 0, . . . , 0) and rank(S') = r. Let M = USV = 
Mq + Mip + . . . + Msp"", where s = C(logp r) and Mi e Zp^". We conjecture 
that when p = 2, 



carry 



A, 



UiSoVo + UoSoVi + UoSiVo + UqSoVo quop 




and when p = 2k + 1, we conjecture that 



ran 




- 2r. 



* Details are outside the scope of presenting this conjecture. 
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Furthermore, in the generic case where U, V are uniformly chosen at random 
in Zp, and n is arbitrarily large, the ranks are equal to bounds above. 

This conjecture shows that a large dimension product of matrices with 
entries in Zp of small rank can still have very large, but not full, rank carry 
matrices. These carries will impact many digits in the expanded product. 
The evidence for the conjecture is experimental^. We developed the formulas 
above by reverse engineering sequences of computed ranks resulting from 
experiments with different primes and matrices of different dimensions. 

In the generic case of equality, the conjecture could be used to generate 
an algorithm for small primes, e.g. when p = 2. However, without further 
refinements, this approach yields an exponential time algorithm which is 
prohibitive. 

4 Conclusion 

We have given efficient algorithms for computing Smith Normal Form over 
two local rings: F[z]/(/'^) and Z/p'^Z. These are useful components for SNF 
algorithms over F[z] and Z, respectively. These algorithms are efficient in 
the black-box model, which means that they are well suited to sparse and 
structured matrices. The integer algorithm is output sensitive. Its memory 
and time usage grows in proportion to the number of nontrivial invariant 
factors. A memory efficient algorithm without that restriction has not been 
found. In addition, we gave a conjecture about the rank of carries resulting 
from multiplying single-digit p-adic matrices. 
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